# ------------------------------------- #
# Replication code for:
#
# Pomeroy, Caleb "The Damocles Delusion: The Sense of Power Inflates Threat  
# Perception in World Politics," International Organization.
#
# This script reproduces the Russian elite survey re-analysis.
# ------------------------------------- #


# --- set your working directory --- #
setwd("~/Downloads/damocles_delusion_replication/")

# --- libraries --- #
library(ggplot2)
library(plyr)
library(texreg)
library(haven)
library(psych)

# --- load data --- #
russia_df <- haven::read_sav("data/russia_elite_data.sav")
dim(russia_df) # the analyses reported in the paper rely on the 2020 wave of the survey (N = 245), using 14 variables 
# you can access the full dataset at https://doi.org/10.3886/ICPSR03724.v8
# here, I provide a subset of the data needed to reproduce the paper's results

# --- prep data for analysis --- #
# --- limited versus extensive Russian national interests
russia_df$extensive_interests <- ifelse(russia_df$FPNATINT == 2, "extensive",
                                        ifelse(russia_df$FPNATINT == 1, "limited", NA))
russia_df$extensive_interests <- factor(russia_df$extensive_interests, levels = c("limited", "extensive"))

# --- militant orientation
russia_df$mil_force_arbiter <- ifelse(russia_df$MILROLE == 1, 1,
                                      ifelse(russia_df$MILROLE == 2, 0, NA))
# --- Ukraine hostile to Russia
russia_df$ukraine_threat <- ifelse(russia_df$UKRAINEF == 98, NA, russia_df$UKRAINEF)

# --- NATO expansion threatens Russian security 
russia_df$nato_threats <- ifelse(russia_df$NATOEXP %in% c(98,99), NA, russia_df$NATOEXP)

# --- US threatens Russian security
russia_df$us_threat_to_natsec <- ifelse(russia_df$FPUSFEAR == 1, 1, 
                                        ifelse(russia_df$FPUSFEAR == 2, 0, NA))
# --- foreign policy elite?
russia_df$elite_type <- factor(russia_df$WORKGROUP)
russia_df$elite_type <- plyr::revalue(russia_df$elite_type, c("1" = "media",
                                                              "2" = "science",
                                                              "3" = "business",
                                                              "4" = "soe",
                                                              "5" = "exec_branch",
                                                              "6" = "leg_fp",
                                                              "9" = "mil_sec"))
russia_df$fp_elite <- ifelse(russia_df$elite_type %in% c("exec_branch", "leg_fp", "mil_sec"), "fp_elite", "other_elite")
russia_df$fp_elite <- factor(russia_df$fp_elite, levels = c("other_elite", "fp_elite"))
table(russia_df$fp_elite) # 105 foreign policy elites, 140 other elites, as footnoted in the main text

# --- gender
russia_df$gender <- ifelse(russia_df$SEX == 1, "male", "female")
russia_df$gender <- factor(russia_df$gender, levels = c("female", "male"))

# --- age
russia_df$age <- (2020-ifelse(russia_df$DOB %in% c(98,99), NA, russia_df$DOB))

# --- prior military service? 
# coded such that higher values = more intensive service (i.e. 3 = served in armed forces, 2 = only in the reserves, 1 = didn't serve at all)
russia_df$mil_service <- (4-ifelse(russia_df$ARMY %in% c(98,99), NA, russia_df$ARMY))

# --- member of United Russia, Putin's de facto party? 
# serves as proxy for support for Putin
russia_df$pid <- ifelse(russia_df$PARTY %in% c(98,99) | russia_df$PARTYMEM %in% c(98,99), NA,
                        ifelse(russia_df$PARTYMEM==1 & russia_df$PARTY==1, "united_party", "other"))
russia_df$pid <- factor(russia_df$pid, levels = c("other", "united_party"))

# --- impact of Russia's foreign policy in recent years on Russia's international influence
# reverse-coded so higher values = perceptions of greater influence
russia_df$russian_influence <- (6-ifelse(russia_df$FPRUINFL %in% c(98,99), NA, russia_df$FPRUINFL))

# --- change in Russia's international influence over last twenty years under Putin
# coded so higher values = perceptions of greater influence
russia_df$russian_influence_putin <- ifelse(russia_df$EVPUINFL %in% c(98,99), NA,
                                            ifelse(russia_df$EVPUINFL == 1, 3,
                                                   ifelse(russia_df$EVPUINFL == 2, 1,
                                                          ifelse(russia_df$EVPUINFL == 3, 2, NA))))

# --- change in Russia's military readiness and strength over last twenty years under Putin
# coded so higher values = perceptions of greater strength
russia_df$military_strength_putin <- ifelse(russia_df$EVPUMIL %in% c(98,99), NA,
                                            ifelse(russia_df$EVPUMIL == 1, 3,
                                                   ifelse(russia_df$EVPUMIL == 2, 1,
                                                          ifelse(russia_df$EVPUMIL == 3, 2, NA))))
# --- reduce power items to a single sense of power scale
# these items serve as a useful proxy for the sense of strength and influence at the core of this paper's concept of the sense of power
ltm::cronbach.alpha(na.omit(cbind(russia_df$russian_influence,
                                  russia_df$russian_influence_putin,
                                  russia_df$military_strength_putin))) 
# cronbach's alpha = 0.70 for the three power items, as reported in the main text
# reduce to unidimensional factor scores
power_bind <- 
  cbind(russia_df$russian_influence,
        russia_df$russian_influence_putin,
        russia_df$military_strength_putin) 
fa.parallel(power_bind)
power_fit <- fa(power_bind, nfactors=1, rotate="varimax", scores=TRUE, fm="pa")
russia_df$sense_power_factor <- power_fit$scores[,1]

# --- regression models --- #
summary(us_main <- glm(us_threat_to_natsec ~ sense_power_factor +
                         extensive_interests + mil_service + mil_force_arbiter + fp_elite + pid + gender + age,
                       data = russia_df, family = "binomial"))
summary(nato_main <- lm(nato_threats ~ sense_power_factor +
                          extensive_interests + mil_service + mil_force_arbiter + fp_elite + pid  + gender + age,
                        data = russia_df))
summary(ukraine_main <- lm(ukraine_threat ~ sense_power_factor +
                             extensive_interests + mil_service + mil_force_arbiter + fp_elite + pid  + gender + age, 
                           data = russia_df))

# --- appendix table A11
screenreg(list(us_main, nato_main, ukraine_main))
# texreg(list(us_main, nato_main, ukraine_main),
#        booktabs = T,
#        caption.above = T,
#        caption = "Sense of Russian Power Explains Threat Perception",
#        custom.coef.names = c("(Intercept)",
#                              "Sense of Russian Power",
#                              "Extensive Russian Interests",
#                              "Prior Military Service",
#                              "Militant Orientation",
#                              "Foreign Policy Elite",
#                              "United Party Member",
#                              "Male",
#                              "Age"),
#        stars = c(.001, .01, .05, .10),
#        symbol = "\\wedge",
#        custom.model.names = c("US Threatens Russian Security",
#                               "NATO Expansion Dangerous to Russia",
#                               "Ukraine Hostile Towards Russia"))

# --- format for plotting purposes
mod_list <- list("us_main" = us_main, 
                 "nato_main" = nato_main,
                 "ukraine_main" = ukraine_main)
results_df <- data.frame()
for(i in 1:length(mod_list)){
  
  mod_i <- mod_list[[i]]
  
  results_df <- rbind(results_df,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
}
results_df$variable <- revalue(results_df$variable, c("sense_power_factor" = "Sense of\nRussian Power",
                                                      "extensive_interestsextensive" = "Extensive\nRussian Interests",
                                                      "mil_service" = "Prior Military\nService",
                                                      "mil_force_arbiter" = "Militant\nOrientation",
                                                      "fp_elitefp_elite" = "Foreign\nPolicy Elite",
                                                      "pidunited_party" = "United Party\nMember",
                                                      "gendermale" = "Male",
                                                      "age" = "Age")) 
levels(results_df$model) <- c("us_main", "nato_main", "ukraine_main")
results_df <- results_df[!results_df$variable%in%c("(Intercept)"),]
results_df$plot_col <- ifelse(results_df$variable == "Sense of\nRussian Power", "#54545d", "#99998e") 
plot_col <- ifelse(results_df$variable == "Sense of\nRussian Power", "#54545d", "#99998e")
results_df$model <- factor(results_df$model, levels = c("us_main", "nato_main", "ukraine_main"))
levels(results_df$model) <- c("US Threatens\nRussian Security", "NATO Expansion\nDangerous to Russia",
                              "Ukraine Hostile\nTowards Russia")
results_df <- results_df[!(results_df$variable %in% c("Age", "Male")),]
results_df$variable <- factor(results_df$variable, levels = rev(c("Sense of\nRussian Power", 
                                                                  "Extensive\nRussian Interests",
                                                                  "Prior Military\nService", 
                                                                  "Foreign\nPolicy Elite", 
                                                                  "United Party\nMember",
                                                                  "Militant\nOrientation")))
# --- figure 3 from the main text
ggplot(results_df, aes(x = est, y = variable, color = plot_col)) +
  geom_vline(xintercept = 0, linetype = "solid", color = "black", size = .8) +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1, color = "black") +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.8) +
  geom_point(size = 5, color = "black") +
  geom_point(size = 4) +
  facet_wrap(~model, nrow = 1, scales = "free_x") +
  theme_bw() +
  labs(x = "Estimate", y = NULL) +
  scale_color_manual(values = plot_col) + 
  theme(axis.title.x = element_text(size = 17, margin = margin(t=15)),
        axis.text.y = element_text(size = 15, color="black"),
        axis.text.x = element_text(size = 13),
        strip.text = element_text(size = 16, margin = margin(r=8,l=8,t=8,b=8)),
        legend.position = "none",
        panel.spacing = unit(1.2, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_rect(fill = "gray90", color = "black", size = 1.5),
        plot.margin = margin(r=75, l=55, t= 15))

# --- demographic distributions for appendix
{plot_demos <- russia_df
  plot_demos$ids <- rownames(plot_demos)
  plot_demos$age_plot <- ifelse(plot_demos$age < 30, "18-29",
                                ifelse(plot_demos$age >= 30 & plot_demos$age < 40, "30-39",
                                       ifelse(plot_demos$age >= 40 & plot_demos$age < 50, "40-49",
                                              ifelse(plot_demos$age >= 50 & plot_demos$age < 60, "50-59",
                                                     ifelse(plot_demos$age >= 60 & plot_demos$age < 70, "60-69", "70+")))))
  plot_demos <- reshape2::melt(plot_demos[,c("ids", "gender", "age_plot", "pid", "fp_elite")], id.vars = "ids")
  plot_demos$variable <- plyr::revalue(plot_demos$variable,
                                       c("gender" = "Gender",
                                         "age_plot" = "Age",
                                         "pid" = "Party ID",
                                         "fp_elite" = "Foreign Policy Elite"))
  plot_demos$ids <- rownames(plot_demos)
  p1 <- 
    ggplot(subset(plot_demos, variable == "Gender"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Female", "Male")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p2 <- 
    ggplot(subset(plot_demos, variable == "Age"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p3 <- 
    ggplot(subset(plot_demos, variable == "Party ID"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Other\nParty", "United\nParty", "NA")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  p4 <- 
    ggplot(subset(plot_demos, variable == "Foreign Policy Elite"), aes(x = value)) +
    geom_bar(color = "black", fill = "gray50", size = .7) +
    facet_wrap(~variable, ncol = 1, scales = "free_x") +
    scale_x_discrete(labels = c("Foreign Policy Elite", "Other Elite")) +
    theme_bw() +
    labs(x = NULL, y = "Count") +
    theme(axis.text.x = element_text(size = 13),
          axis.text.y = element_text(size = 12),
          axis.title.y = element_text(size = 13),
          strip.text.x = element_text(size = 13))
  }

# --- appendix figure A3
gridExtra::grid.arrange(p1,p2,p3,p4, ncol = 1)

